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I.  INTRODUCTION 


A.  BACKGROUND 

A  primary  goal  of  a  computer  graphics  system  is  to  provide  the  user  with 
different  views  of  objects.  Sometimes  the  objects  that  the  user  manipulates  are 
simple  in  nature  and  can  be  constructed  easily  with  the  primitives  provided  by 
the  graphics  support  package.  However,  in  most  applications  areas  the  user  is 
concerned  with  more  complex  objects.  The  display  of  three-dimensional  surfaces  is 
one  such  application  area  in  computer  graphics  and  is  the  area  that  this  study 
explores. 

It  is  the  primary  intent  of  this  study  to  stimulate  the  reader’s  interest  in  the 
area  of  three-dimensional  surface  generation  and  display.  To  provide  this 
stimulation,  we  combine  the  power  of  certain  mathematical  techniques  and  a  high 
performance  graphics  environment  to  design  and  implement  a  set  of  functions 
that  can  be  used  to  create,  manipulate,  and  display  three-dimensional  solid-filled 
surfaces.  Once  developed,  the  reader  will  not  only  be  able  to  use  these  functions 
to  explore  the  design,  representation,  and  rendering  of  such  surfaces  but  also  will 
be  able  to  use  these  functions  in  other  fields  that  can  derive  benefit  from  their  use 
such  as  cartography,  robotics,  computer  vision,  and  artificial  intelligence. 
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B.  PROJECT  DESCRIPTION 


1.  Motivation 

For  developing  graphics  applications,  the  typical  graphics  environment 
provides  the  user  with  a  variety  of  sophisticated,  powerful  tools  and  a  high-level 
language.  Almost  inevitably,  an  application  calls  for  something  that  is  not 
provided  for  or  if  provided  is  not  acceptable.  The  lack  of  solid-filled  curved 
surface  support  in  one  particular  high-performance  graphics  workstation 
prompted  this  study.  By  providing  this  capability  we  can  create  applications  that 
enhance  and  expand  what  we  are  able  to  do  and  conceive  of  doing  on  these 
workstations. 

2.  Proposed  Capabilities 

To  provide  this  capability  for  generating  solid-filled  parametric  bicubic 
surface  patches,  we  must  look  for  a  method  that  is  simple,  powerful, 
understandable,  and  implementable  in  software.  One  approach  would  be  to  start 
from  scratch  and  develop  a  solution.  This,  however,  is  usually  not  a  wise  way  to 
proceed.  A  solution  derived  in  this  manner  is  likely  to  have  severe  limitations 
such  as  being  computationally  inefficient,  lacking  robustness,  and  difficult  to  use. 
Another  more  prudent  way  to  proceed  is  to  find  some  basic  mathematical 
techniques  that  can  be  applied  to  the  problem  in  a  form  consistent  with  the 
existing  graphics  environment.  Solutions  developed  in  this  manner  are  more 
easily  accepted  since  their  basis  lies  in  proven  mathematical  techniques  and 
generally  meets  the  requirements  we  seek  of  being  simple,  powerful, 
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understandable  and  implementable.  We  chose  to  develop  our  new  capability 
following  the  latter  method  -  using  pre-existing  mathematics. 

a.  Parametric  Bicubic  Surface  Construction 

One  of  the  primary  capabilities  needed  when  working  with  surfaces  is 
to  have  an  effective  method  for  describing,  manipulating  and  displaying  them. 
For  this,  we  base  our  work  with  surfaces  on  a  mathematical  model  -  the 
parametric  bicubic  surface  patch.  Choice  of  this  method  allows  the  user  to  rapidly 
develop  a  smooth  surface  and  display  it  as  a  wireframe  image  through  the 
specification  of  only  a  few  points  called  control  points. 

b.  Polygonal  Parametric  Bicubic  Surface  Decomposition 

Another  capability  required  is  to  extract,  from  the  mathematical 
model,  the  information  needed  to  construct  a  solid  surface.  For  this,  we  need  to  be 
able  to  decompose  any  surface  represented  in  parametric  bicubic  form  into 
arbitrarily  small  polygons.  These  polygons  can  then  be  used  in  the  construction 
of  a  polygon  mesh.  This  capability  permits  the  user  to  use  standard  or 
customized  algorithms  to  manipulate  the  surface. 

C.  PROGRAMMING  ENVIRONMENT 

The  IRIS  Turbo  2400  Graphics  Workstation,  manufactured  by  Silicon 
Graphics,  Inc.,  is  the  target  programming  environment  for  the  design, 
development,  and  implementation  of  the  parametric  bicubic  surface  constructor 
and  polygonal  bicubic  surface  decomposition  functions.  The  IRIS  has  special- 
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purpose  hardware  that  is  designed  to  replace  less  efficient  software.  The  system 
supports  real-time  color  graphics,  the  Unix  operating  system  software,  and  the  C 
programming  language.  A  high  resolution  color  monitor  provides  the  output 
device  for  displaying  the  graphical  output  of  the  modeling  and  user  defined 
functions. 

In  addition  to  the  standard  programming  support  provided  by  the  UNIX 
operating  system,  the  IRIS  has  an  extensive  collection  of  utility  and  graphics 
functions  contained  in  a  Graphics  Library.  This  library  provides  the  user  high- 
level  access  to  the  hardware,  enabling  graphical  objects  to  be  easily  manipulated 
as  geometrical  objects  (points,  lines,  polygons,  etc.)  rather  than  pixels.  A  series  of 
coordinate  systems  and  mapping  instructions  also  provides  the  user  with  the 
capability  to  define  such  objects  in  world  coordinate  space. 


10 


II.  THREE  DIMENSIONAL  REPRESENTATION  OF  SURFACES 


We  stated  in  the  previous  chapter  that  an  important  application  area  in 
computer  graphics  is  concerned  with  the  three-dimensional  representation  of 
surfaces  but  we  did  not  describe  how  one  might  do  this.  There  are  many  ways  to 
represent  surfaces.  The  two  most  commonly  used  representations  are:  polygon 
meshes  and  parametric  bicubic  patches. 

A.  POLYGON  MESHES 

A  polygon  mesh  is  nothing  more  than  a  set  of  connected  polygonally  planar 
surfaces.  There  are  several  ways  in  which  these  planar  polygons  can  be 
represented  -  explicit  polygons,  vertex  list  pointers,  explicit  edges,  etc..  Each  of 
these  representations  has  its  own  advantages  and  disadvantages  and  various 
criteria  can  be  applied  to  evaluate  the  representation.  Criteria  such  as  how  much 
primary  and  secondary  storage  is  available,  how  easy  is  it  to  identify  the  polygons 
sharing  an  edge,  how  difficult  is  it  to  display  the  mesh,  name  only  a  few 
commonly  used  for  evaluation.  Regardless  of  how  the  mesh  is  represented,  there 
are  many  algorithms  available  for  processing  it.  For  example,  algorithms  have 
been  developed  to  remove  hidden  surfaces  from  an  object,  while  others  can 
produce  lighting  and  shading  effects  on  the  surface  of  the  object. 
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Many  objects  that  have  planar  features  such  as  buildings,  tables,  and  cabinets 
can  easily  be  represented  by  a  polygon  mesh.  However,  polygon  meshes  are  not 
limited  to  representing  only  planar  objects.  They  can  also  be  used  to  represent 
curved  surfaces.  To  represent  a  curved  surface  with  a  polygon  mesh,  one  creates 
an  approximation  to  the  surface  by  using  arbitrarily  small  polygons.  The  smaller 
the  polygons  the  better  the  approximation.  However,  certain  difficulties  arise 
when  representing  surfaces  as  polygon  meshes.  As  one  approximates  the  surface 
with  smaller  and  smaller  polygons,  both  space  and  execution  time  of  the 
algorithms  that  process  the  mesh  increase  linearly. 

B.  PARAMETRIC  BICUBIC  PATCHES 

Parametric  bicubic  patches  are  mathematical  models  of  a  surface  and 
represent  one  of  the  simplest  mathematical  elements  we  can  use  to  model  an 
arbitrary  surface.  A  patch  can  be  defined  as  a  curve  bounded  collection  of  points 
whose  coordinates  are  given  by  a  continuous,  two-parameter,  single  valued 
function  in  the  form 

x  =  x(s,t) 

y  =  yM) 

z  =  z(s,t) 

where  the  parameters  are  restricted  to 

s,t  e  [0,1]. 
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Patches  are  used  to  describe  the  surface  of  an  object  in  a  piecewise  manner 
much  like  the  polygon  mesh.  That  is,  many  patches  are  used  to  describe  the 
surface  of  an  object  and  the  total  surface  is  generated  by  displaying  all  of  its 
individual  patches. 

Bicubic  patches  are  most  often  used  to  generate  line  drawings  of  three- 
dimensional  objects  -  often  called  wireframe  representations.  One  benefit  of  using 
bicubic  patches  is  that  fewer  bicubic  patches  than  polygonal  patches  are  needed 
to  represent  a  curved  surface  to  a  given  accuracy.  On  the  other  hand,  the 
algorithms  for  working  with  bicubic  surfaces  are  more  complex  than  those  for 
polygon  meshes.  In  addition,  wireframe  representations  can  exhibit  certain 
deficiencies.  For  example,  objects  that  are  represented  by  wireframes  are  often 
ambiguous.  That  is,  you  can  look  at  the  same  object  and  get  different  visual 
interpretations.  One  reason  for  the  different  interpretations  is  the  ability  to  see 
through  a  wireframe  representation  of  an  object.  While  not  a  problem  for  some 
applications,  it  can  be  a  serious  problem  for  others. 
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III.  PARAMETRIC  CUBIC  CURVES  AND  SURFACES 


A.  GENERAL 

The  method  that  we  use  to  generate  surfaces  is  based  on  using  parametric 
cubic  curves  so  it  is  helpful  to  review  the  mathematical  basis  for  these  curves.  A 
parametric  cubic  curve  has  the  property  that  x,  y,  and  z  can  be  defined  els  third- 
order  polynomials  for  some  variable  t  [Ref.  l:pp.  34]  : 

x(t)  =  axts  +  bxt2  +  cxt  +  dx, 
y(t)  =  ayts  +  byt2  +  cyt  +  dy, 
z(t)  =  atts  +  bxt2  +  ctt  +  dt. 

These  equations  are  known  as  the  algebraic  form  of  a  parametric  cubic  curve. 
In  this  form,  we  can  identify  a  unique  set  of  12  constant  algebraic  coefficients. 
These  algebraic  coefficients  determine  a  unique  parametric  cubic  curve;  they 
determine  the  size  and  shape  of  the  curve  and  its  position  in  three-dimensional 
space.  Two  curves  of  the  same  shape  have  different  algebraic  coefficients  if  they 
occupy  different  positions  in  three-dimensional  space.  Because  we  want  to  deal 
with  a  finite  segments  of  the  curve,  we  limit  the  range  of  the  parameter,  without 
loss  of  generality,  to  O^t^l.  We  call  these  finite  pieces  curve  segments.  A  curve 
segment  is  nothing  more  than  a  point-bounded  collection  of  points.  In  our  case 
the  points  are  bounded  at  t=0  and  t=l.  Each  equation  in  the  algebraic  form  can 
be  expressed  as  a  vector  product  as  follows: 
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x(t)  =  [t*  t2  t  l] 


a 


=  [t*t2t  l]cx  =  TCX. 

Using  the  vector  product  form  is  usually  more  convenient  as  it  separates  the 
distinct  parameters  of  the  parametric  equation  into  unknown  coefficients  of  x(t) 
and  the  parameter  t  that  we  wish  to  manipulate.  Here  T  is  the  row  vector  of 
powers  of  t,  while  Cx  is  the  column  vector  of  coefficients  of  x(t).  Similarly  the 
parametric  equations  for  y(t)  and  z(t)  can  be  written  as  y(t)  =  TCy  and  z(t)  =  TCt. 
By  varying  the  parameter  t  from  0  to  1  in  each  equation  we  define  the  curve 
segment. 

Arbitrarily  assigning  values  to  these  unknown  coefficients  results  in  defining  a 
curve  in  three-dimensional  space.  However,  it  is  not  easy  to  determine  the 
properties  of  this  curve.  What  we  wish  to  do  is  establish  some  constraints  on 
these  coefficients.  We  want  the  curves  we  generate  to  have  some  predictable 
properties.  To  solve  the  equations  for  these  unknown  algebraic  coefficients,  we 
establish  a  set  of  constraints,  thereby  defining  a  unique  cubic  curve  with 
predictable  properties.  To  illustrate  this  process  we  look  at  some  example  cubic 
curves  for  which  the  constraints  are  well-known. 

B.  CUBIC  CURVE  EXAMPLES 

1.  Hermite  Curve 

The  Hermite  cubic  curve  is  determined  from  its  endpoints  (Pt,  P2)  and 
endpoint  tangents  (Rlr  R2).  In  the  literature  [Ref.  2:pp.  516-519]  [Ref.  3:pp.  123- 
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129],  we  find  the  geometric  form  of  the  Hermite  curve  to  be 

Qh(t)  =  TMhGh. 

In  geometric  form.  T  is  the  row  vector  of  the  powers  of  the  parametric 
variable  t,  Mh  is  the  basis  matrix,  and  Gh  is  the  geometry  vector.  A  basis  matrix 
refers  to  a  constraint  procedure  that  is  embodied  in  matrix  form  and  a  geometry 
vector  contains  the  control  points  used  to  guide  the  curve. 

In  this  particular  case,  the  Hermite  basis  matrix  (Mh)  is 

2-211 
-3  3  -2  -1 

0  0  10 
10  0  0 

and  the  Hermite  geometry  vector  is 

Pi 

P2 

R,  • 

Rj 

Now  using  the  above  formulation,  given  two  points  and  their  tangents,  we 
can  evaluate  x(t),  y(t),  and  z(t)  for  0^t$l  and  find  all  points  on  the  Hermite  form 
of  the  cubic  curve  from  Pi  to  P2  with  starting  tangent  vector  R,  and  ending 
tangent  vector  R2.  It  is  through  these  constraints  (Mh)  that  the  control  points 
(Gh)  control  the  parametric  equations  and  produce  an  equation  that  can  generate 
a  discretely  sampled  curve  segment  in  three-dimensional  space. 
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If  we  take  the  product  TMh,  we  have 


TMh  =  [(2ts-3t2  +  l)  (  — 2ts  +  3t2)  (ts-2t2  +  t)  (t*-t2)  j. 

These  four  functions  of  t  in  the  product  TMh  are  often  called  blending  functions 
[Ref  l:pp.  48-52].  As  the  name  implies,  they  blend  the  effects  or  contributions  of 
the  endpoints  and  tangent  vectors  to  produce  the  intermediate  point  coordinate 
values  over  the  domain  of  t. 

2.  Bezier  Curve 

The  defining  form  for  a  Bezier  cubic  curve  is  similar  to  the  Hermite  form. 
The  difference  is  in  the  definition  of  the  endpoint  tangent  vectors.  The  Bezier 
form  uses  four  points  (  P19  P2.  Ps,  P4  )  instead  of  2  points  and  2  tangent  vectors. 
The  tangent  vectors  at  the  endpoints  in  Bezier  form  are  determined  by  the  line 
segments  P^  and  PSP4  .  The  Bezier  cubic  curve  passes  through  the  first  and 
fourth  control  points  (P,  and  P4)  and  uses  the  second  and  third  points  (P2  and 
Ps)  to  determine  the  shape  of  the  curve.  [Ref.  l:pp.  113-125]  [Ref.  2:pp.  519-521]. 
The  geometric  form  of  the  Bezier  curve  is 


Qb(t)  =  TMbGb 


where  the  Bezier  basis  matrix  (Mb)  is 


Mb  — 


-1  3  -3 

3-6  3 

-3  3  0 

10  0 


1 

0 

0 

0 
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and  the  Bezier  geometry  vector  is 


Pi 

P2 

Ps  • 

Pi 

The  Bezier  form  of  the  cubic  curve  is  more  widely  used  in  computer 
graphics  than  the  Hermite  form.  A  primary  reason  for  its  popularity  is  that  the 
geometry  matrix  of  four  points  (Gb)  is  more  intuitive  for  an  interactive  user.  The 
user  has  only  to  manipulate  the  four  points  and  does  not  have  to  specify  the 
tangent  vectors.  It  is  usually  easier  for  a  person  to  think  about  manipulating 
points  rather  than  trying  to  manipulate  points  and  tangent  vectors. 

3.  Other  Useful  Cubic  Curves 

The  Hermite  and  Bezier  cubic  curves  are  not  the  only  forms  of  cubic 
curve  that  are  available.  Two  others  are  the  Cardinal  Spline  and  the  B-Spline. 

a.  Cardinal  Spline 

The  Cardinal  Spline  curve  passes  through  the  two  interior  control 
points  (P2  and  P8)  and  uses  the  points  Pj  and  P4  to  define  the  shape  of  the  curve 
[Ref.  4:p.  11-4].  The  geometric  form  of  the  Cardinal  curve  is 

Qe(t)  =  TMeGc 

where  the  Cardinal  basis  matrix  (Mc)  is 
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—  a  2  — a  a  — 2  a 

2a  a  —  3  3  2a  a 

—a  0  a  0 

0  10  0 

and  the  Cardinal  geometry  vector  is 

Pi 

P2 

P* 

P4 

The  scalar  coefficient  a  in  the  Cardinal  basis  matrix  must  be  positive  and 
determines  the  length  of  the  tangent  vector  at  point  P2  and  Ps. 
b.  B-Spline 

The  geometric  form  of  the  B-Spline  curve  is 

Qb«(t)  =  TMb,Gb8 

where  the  B-Spline  basis  matrix  (Mbs)  is 

-1  3-31 

3-630 
-3  0  3  0 

14  10 

and  the  B-Spline  geometry  vector  is 

Px 

P2 

P*  • 
p< 
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In  general,  the  B-Spline  curve  does  not  pass  through  any  control 
points  but  is  continuous  and  also  has  continuity  of  tangent  vectors  and  of 
curvature  (that  is,  first  and  second  derivatives  are  continuous  at  the  endpoints). 
The  Hermite  and  Bezier  forms  have  only  first-derivative  continuity  at  the 
endpoints,  but  do  pass  through  control  points  [Ref.  l:pp.  125-146]  [Ref.  2:pp. 
521-523]. 

C.  DEFINING  SURFACES 

By  adding  a  new  parameter  s  and  additional  algebraic  coefficients  to  the  cubic 
curves  in  the  previous  section,  we  can  define  the  algebraic  form  of  a  bicubic 
surface  patch  [Ref.  l:pp.  156]  as 

i=s  j=s 

P(s>t)  =  E  Eaus'tJ 

1=0  j=o 

with  the  restriction  on  the  parametric  variables  to 

S,t  £  [0,1]. 

By  varying  both  parameters  from  0  to  1  in  each  equation,  we  define  all 
points  on  the  surface  patch.  Assigning  one  parameter  a  constant  value  and 
varying  the  other,  results  in  a  cubic  curve. 

Expanding  the  above  equation  in  terms  of  x(s,t)  and  noting  that  the  terms  for 
y(s,t)  and  z(s,t)  are  similar  we  have 
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x(s.t)  =  a88s8t8  +  a82sst2  +  a81s*t  +  a80ss  +  a28s2t8  +  a22s2t2  +  a21s2t  +  a20s2 
+  aiss^*  +  ai2st2  +  anst  +  aios  +  a0  jt*  +  a02^2  +  aOlt  +  a00" 

Written  in  vector  product  form 


x(s,t)  =  SC„TT 


where 


S  =  |ss  s2  s  1  J , 
T  =  [t8  t2  t  1  ]  , 


aS2 

aSl 

aso 

1 

a2S 

a22 

a21 

a20 

X  “ 

al8 

a  12 

all 

a10 

aos 

a02 

a01 

o 

o 

CQ 

and  TT  is  the  transpose  of  the  matrix  T  . 

From  these  equations  we  can  see  that  there  are  48  degrees  of  freedom  or 
algebraic  coefficients  that  we  must  specify.  Like  the  cubic  curve,  a  change  in  any 
one  of  these  coefficients  defines  a  different  surface. 

The  complete  algebraic  manipulation  of  the  equations  to  arrive  at  the 
following  equation  is  similar  to  that  of  the  curve  process  described  in  the  previous 
section.  For  the  Bezier  surface  patch,  the  geometric  form  of  the  equation  is: 

x(s,t)  =  SMbQxMjTT 

where  Mb  is  the  same  Mb  as  in  the  Bezier  curve  equation,  Mb  is  its  transpose,  and 
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Qx  is  the  x  component  of  the  sixteen  control  points  of  a  surface  patch.  The 
matrix  Qx  is 


Qx  = 


and  similarly  for  Qy  and  Qr. 

Since  we  must  provide  three  4x4  matrices,  one  for  each  of  component  x,  y, 
and  z,  it  can  be  seen  that  we  have  specified  the  48  degrees  of  freedom  as  in  the 
algebraic  form. 

As  can  be  seen  by  the  above  equations,  a  bicubic  surface  patch  can  be  defined 
by  a  set  of  16  control  points  and  a  basis  matrix.  By  manipulating  the  control 
points,  we  change  the  shape  of  the  surface  as  constrained  by  the  basis  matrix. 
We  take  this  knowledge  with  us  as  we  design  our  functions  in  the  next  chapter. 
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IV.  DESIGN  AND  IMPLEMENTATION  OF  OUR  SURFACE  FUNCTIONS 


A.  OVERVIEW 

Having  the  mathematics  developed  in  the  previous  chapter  does  us  no  good  if 
we  can  not  to  put  it  to  use.  In  our  case,  this  means  being  able  to  draw  a 
parametric  bicubic  surface.  It  is  at  this  point  that  we  begin  to  see  how  the 
mathematics  can  be  combined  with  the  power  of  the  computer  and  the  graphics 
workstation. 

Up  to  this  point,  we  have  dealt  with  two  forms  of  the  parametric  cubic  curve 
and  parametric  bicubic  surface  -  the  algebraic  form  and  the  geometric  form.  The 
question  now  is  which  one  shall  we  work  with? 

Deciding  what  form  to  use  depends  largely  on  the  application.  If  we  are  given 
or  know  the  algebraic  equations  of  the  curve,  then  the  reasonable  choice  is  the 
algebraic  form.  If  we  plan  to  do  surface  fitting  of  data  or  interactive  design  the 
choice  is  the  geometric.  We  choose  to  use  the  geometric  form.  Our  primary' 
reason  for  choosing  it  is  that  the  geometric  form  offers  us  a  greater  insight  into 
the  control  and  behavior  of  curves  and  surfaces  than  is  otherwise  available  with 
the  classical  algebraic  formulation.  It  should  be  noted,  however,  that  it  is 
possible,  through  mathematical  manipulation,  to  convert  from  one  form  to  the 
other  [Ref.  l:pp.  164]. 
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B.  METHODOLOGY 


How  one  proceeds  to  generate  a  surface  impacts  usefulness,  flexibility,  and  the 
ability  to  understand.  As  we  stated  in  the  first  chapter,  the  built-in  functions 
provided  by  a  graphics  environment  are  not  always  exactly  what  we  want.  This 
is  the  case  in  the  IRIS  graphics  environment.  Although  the  IRIS  provides  support 
for  bicubic  surface  patches,  it  does  not  support  surface  patch  decomposition. 
What  we  want  is  the  capability  to  do  both.  We  also  want  this  capability  without 
sacrificing  what  a  user  already  knows  about  how  the  IRIS  supports  bicubic  surface 
patches.  To  achieve  this,  we  have  developed  a  set  of  parallel  routines  [Appendix 
A]  that  provide  nearly  all  the  functionality  as  the  standard  IRIS  functions  while 
at  the  same  time  providing  the  user  with  extended  support  via  three  additional 
functions.  These  extended  functions  allow  the  user  to  have  access  to  the 
triangular  polygons  that  our  new  functions  generate  during  the  construction  of 
the  bicubic  surface  patch.  That  is,  the  IRIS  user  is  able  to  use  the  new  functions 
in  the  same  way  as  he  would  use  the  standard  functions  by  substituting  the 
names  of  the  new  routines  in  place  of  the  standard  IRIS  routines.  If,  however,  the 
user  wishes  to  be  able  to  have  access  to  the  individual  triangular  polygons  that 
make  up  the  surface,  he  has  only  three  additional  routines  to  learn. 

C.  PARALLELING  THE  IRIS  SUPPORTED  FUNCTIONS 

The  IRIS  graphics  environment  has  available  five  functions  for  defining  and 
generating  parametric  bicubic  surface  patches.  Those  five  functions  are: 
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-  dcfl>asis  defines  a  basis  matrix 

-  patchbasis  sets  the  current  basis  matrices  for  both  the  s  and  the  t  parametric 
direction 

-  patchcurvcs  sets  the  number  of  curves  used  to  represent  a  patch 

-  patchprccision  sets  the  precision  at  which  the  curves  are  drawn 

-  patch  draws  the  surface  patch. 

A  complete  description  of  the  functions  and  their  arguments  can  be  found  in 
the  IRIS  Users  Manual  [Ref.  4]. 

To  ease  the  pain  of  learning  new  functions,  our  new  functions  are 
syntactically  identical  to  the  standard  IRIS  functions  with  the  exception  that  the 
new  function  names  are  the  standard  function  names  prefixed  by  the  letter  n. 
These  parallel  functions  are: 


-  ndcfbasis 

-  npatchbasis 

-  npatchcurvcs 

-  npatchprccision 

-  npatch . 

The  usage  and  the  arguments  of  these  parallel  functions  remain  the  same  as 
the  standard  IRIS  functions.  The  only  difference  that  the  user  notices  is  that  the 
wireframe  drawn  looks  like  a  triangular  mesh  instead  of  the  typical  wireframe  and 
that  the  function  npatchprccision  has  no  effect  on  the  displayed  image. 

While  these  routines  seem  to  do  what  the  old  routines  do,  they  are  more 
powerful  because  they  provide  special  extensions  to  the  user.  These  extensions 
provide  the  capability  to  manipulate  the  surface  patch  as  individual  polygons. 
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D.  TRIANGULAR  DECOMPOSITION  OF  SURFACE  PATCHES 


The  extensions  mentioned  in  the  previous  section  are  available  to  the  user  by 
using  three  additional  routines: 

-  Set  User  Routine  for  npateh  provides  an  intercept  function  for  handling  the 
triangular  polygons  generated  in  the  surface  decomposition 

-  User  Routine  provides  the  user  a  way  to  turn  the  intercept  function  on  and 
off 

-  Set  Default  Routine  for  npateh  allows  the  user  to  return  to  the  system 
defined  intercept  function. 

There  is  one  argument  to  the  function  Set  User  ^Routine  for  npateh.  This 
argument  is  the  name  of  a  user-defined  function  that  expects  to  receive  a  3x3 
array.  This  3x3  array  contains  the  three  vertices  of  a  triangle  where  each  vertex  is 
made  up  of  an  x,  y,  and  z  coordinate.  The  function  User  ^Routine  expects  one 
argument  also.  If  this  argument  is  zero,  then  the  intercept  function  is  turned  off; 
i.e.,  the  user’s  program  cannot  intercept  the  triangles  composing  the  surface. 
Otherwise  the  function  is  activated  allowing  the  user’s  program  to  intercept  the 
triangles  comprising  the  surface.  The  function  Set  Default  Routine  for  npateh 
does  not  expect  any  arguments. 

Using  these  functions,  the  user’s  program  has  access  to  and  can  manipulate 
the  individual  triangular  components  of  the  surface  patch.  For  example,  an 
individual  surface  patch  can  be  decomposed  into  triangular  polygons  and  then  via 
the  user  intercept  function,  each  polygon  can  be  subjected  to  an  illumination 
model  that  produces  a  realistic  looking  surface  in  three-dimensional  space. 
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The  surface  patch  is  decomposed  into  triangular  polygons  one  at  a  time.  For 
the  user's  program  to  intercept  these  polygons,  the  program  must  have  specified 
an  intercept  function  via  the  Set  User  Routine  for  npateh  and  must  have 
activated  it  via  User  Routine.  Then,  as  each  polygon  is  generated,  the  user's 
function  can  process  them  in  any  way  desired.  They  can  be  stored,  manipulated, 
altered,  etc..  It  is  the  user’s  program  that  determines  what  to  do  with  them. 
This  feature  provides  a  tremendous  amount  of  flexibility,  creativity,  and 
applicability  above  what  is  currently  available  in  the  standard  IRIS  graphics 
environment  support  of  surfaces. 

E.  GENERAL  GUIDELINES  FOR  USAGE 

To  prevent  any  unnecessary  problems  in  using  our  new  functions,  we  need  to 
establish  a  basic  set  of  guidelines  or  sequences  of  events  that  should  be  followed. 
If  the  user  does  not  want  to  use  the  special  extensions,  i.e. 
Set  User  Routine  for  npateh  and  User  Routine ,  then  a  modified  version  of  the 
standard  IRIS  setup  steps  for  using  surface  patches  can  be  followed.  These  steps 
are: 

-  define  the  appropriate  curve  bases  using  the  ndefbasis  function; 

-  select  a  basis  for  the  s  and  t  parametric  directions  using  the  npatehbasis 
function; 

-  select  the  number  of  curve  segments  to  be  drawn  in  each  parametric  direction 
using  the  npateheurves  function; 

-  draw  the  surface  by  using  the  npateh  function. 

The  only  change  to  the  standard  IRIS  setup  is  that  it  is  not  necessary  to  use 
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the  npatehpreeision  function.  This  function  does  not  effect  the  displayed  image 
and  its  only  purpose  is  to  maintain  consistency  with  the  standard  IRIS  functions. 

If  the  user  wishes  to  use  the  extensions  via  the  Set  User  Routine  for  npateh 
and  the  User  Routine  functions,  then  the  following  steps  must  added: 

-  An  intercept  function  must  be  declared  and  defined  before  calling  the 
Set  User  Routine  for  npateh  function.  This  intercept  function  must  be 
declared  as  a  function  returning  an  integer  value  (even  though  it  is  not  used) 
and  must  be  defined  as  receiving  a  3x3  matrix  of  floating  point  numbers, 
where  each  row  contains  one  set  of  x,  y  and  z  coordinates  of  an  intercepted 
triangles  vertex.  The  name  of  this  intercept  function  will  be  the  argument 
given  to  the  Set  User  Routine  for  npateh  function; 

-  Activating/deactivating  the  intercept  function  via  the  User  Routine  function 
can  be  performed  any  time  after  the  above  step  has  been  completed. 

By  using  these  guidelines,  a  user  should  not  have  any  difficulty  in  using  these 
functions.  As  we  will  show  in  the  next  chapter,  these  functions  are  easy  to  use, 
efficient,  and  can  provide  some  impressive  results  when  a  carefully  chosen 
intercept  function  is  used. 
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V.  USAGE  AND  PERFORMANCE 


Having  taken  a  brief  tour  of  the  functions  we  designed  in  the  previous 
chapter,  we  need  to  provide  some  concrete  examples  of  their  usage,  performance 
levels,  and  limitations. 

A.  SAMPLE  PROGRAMS 

Appendix  B  gives  the  listing  for  four  sample  programs  using  the  new 
functions.  Each  program  illustrates  how  the  new  functions  can  be  integrated  into 
the  IRIS  graphics  programming  environment.  Program  #1  draws  2  surface 
patches  in  wireframe  representation.  One  surface  is  drawn  using  the  standard 
IRIS  functions  while  the  other  is  drawn  using  the  parallel  functions.  When  this 
program  is  run  the  user  notices  the  different  appearance  of  the  wireframe  surface 
patch  drawn  with  the  new  functions.  It  has  the  triangular  mesh  appearance 
described  in  the  previous  chapter.  Program  #2  shows  how  a  user-defined 
intercept  function  can  be  used  via  the  three  extension  functions  we  have  designed. 
This  program  intercepts  the  triangles  generated  during  the  patches  decomposition 
and  puts  them  into  an  IRIS  graphical  object.  Program  # 3  shows  how  the  user- 
defined  intercept  function  can  be  dynamically  changed  during  execution  and 
Program  #4  shows  how  a  well  chosen  intercept  function  can  be  used  to  produce 
realistic  lighting  of  a  surface. 
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B.  PERFORMANCE  COMPARISONS 


Knowing  that  our  functions  work,  we  would  like  to  know  how  efficient  they 
are.  The  way  that  we  approach  this  question  is  to  compare  our  new  functions  to 
the  standard  IRIS  functions.  Because  we  have  designed  our  functions  to  be 
substitutable  as  a  set  for  the  standard  functions,  we  do  not  have  any  problems  in 
testing  the  relative  performance  of  the  sets.  However,  two  words  of  caution  are  in 
order  before  comparing  these  two  sets  of  functions.  First,  the  IRIS  implements 
many  of  its  graphics  primitives  via  special  purpose  hardware.  This  is  the  case 
with  the  function  patch.  Therefore,  its  parallel  function  npatch,  which  is 
implemented  in  software  is  not  as  fast.  Second,  the  new  function  npatchprccuion 
does  not  affect  the  computation  whereas  the  IRIS  function  patchprccision  does 
affect  the  computations.  Keeping  these  points  in  mind,  we  have  developed  a 
simple  benchmark  program,  listed  in  Appendix  C,  that  we  use  to  draw  a 
wireframe  representation  of  a  surface  patch  100  times.  By  executing  this  program 
10  times  and  getting  the  average  times,  we  can  get  an  idea  of  the  performance  of 
the  parallel  set  of  surface  patch  functions  as  compared  to  the  standard  IRIS 
surface  patch  functions. 

The  way  that  we  measure  performance  of  a  particular  program  is  to  use  the 
UNIX  time  command.  The  time  command  returns,  on  program  completion,  the 
time  in  seconds  for  system  time,  user  time,  and  elapsed  time. 

The  benchmark  program  was  executed  in  two  different  modes.  In  the  first 
mode,  the  program  was  executed  without  the  assistance  of  the  IRIS’s  floating 
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point  accelerator  (FFP)  while  in  the  second  mode,  the  program  was  executed  with 
the  FFP.  The  results  indicate  that  the  standard  IRIS  functions  are  350%  faster 
without  the  FFP  and  260%  faster  with  the  FFP.  As  expected,  the  new  functions 
are  slower,  however  considering  that  they  were  not  designed  to  replace  the  IRIS 
functions,  these  results  are  good.  Normally,  one  can  expect  an  order  of  magnitude 
increase  in  performance  when  special  hardware  is  used. 

C.  LIMITATIONS 

Nothing  that  can  be  developed  is  without  limitations.  The  reader  should 
recall  that  it  was  certain  limitations  of  the  existing  IRIS  system  that  motivated 
this  study.  The  functions  we  have  designed  and  have  implemented  have  allowed 
us  to  overcome  certain  limitations  in  the  IRIS  graphics  environment.  At  the  same 
time,  these  functions  have  their  own  limitations. 

The  primary  limitation  of  our  parallel  functions  is  speed.  While  these 
functions  have  been  carefully  implemented  using  efficient  algorithms  and  data 
structures,  they  are  not  as  fast  as  using  special  purpose  hardware.  Another 
limitation  deals  with  the  use  of  memory.  The  npateh  function  allocates  memory 
to  save  each  point  on  the  surface  patch.  The  number  of  points  that  are  generated 
are  proportional  to  the  product  of  the  desired  number  of  curve  segments  in  the  s 
and  t  parametric  directions.  For  example,  to  draw  a  surface  patch  with  10  curves 
in  the  s  direction  and  10  curves  in  the  t  directions  requires  at  least  enough 
memory  to  store  300  floating  point  numbers  (one  for  each  x,  y,  and  z  component). 
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To  draw  a  100x100  surface  patch  requires  enough  memory  to  store  30,000  floating 
point  numbers.  Assuming  a  floating  point  number  requires  4  bytes,  the  10x10 
patch  requires  1.17  Kilobytes  of  memory,  while  the  100x100  patch  requires  117.1 
Kilobytes  of  memory. 
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VI.  RECOMMENDATIONS  AND  CONCLUSIONS 


A.  DIRECTIONS  FOR  FURTHER  STUDY 

Bicubic  surface  display  and  generation  is  an  area  of  research  in  computer 
graphics  that  is  exciting  and  important.  Since  the  development  of  high- 
performance  graphics  systems,  the  demand  for  realism  and  real-time  has  increased 
significantly.  Consequently,  the  need  is  great  for  continued  creativity  and 
exploration  in  the  area. 

1.  Development  of  Application  Programs 


The  power  of  these  parallel  surface  functions  we  have  created  can  only  be 
derived  through  the  use  of  the  intercept  functions.  Whether  they  will  be  used  to 
experiment  with  lighting  and  shading  models  or  applied  to  fractal  geometry  can 
only  be  answered  by  time.  However,  it  is  through  creative  experimentation  that 
these  questions  can  be  answered.  Some  areas  for  further  work  are: 

-  Surface-fitting  sampled  data 

Surface-fitting  is  the  process  of  constructing  a  representation  to  model  the 
surface  of  an  object  based  on  a  fairly  large  number  of  given  data  points.  By 
taking  these  points  and  chosing  an  appropriate  set  of  surface  constraints,  one 
can  accurately  reconstruct  the  surface.  For  example,  during  this  work,  we 
were  given  a  set  of  digitized  x,  y  and  z  coordinates  for  a  human  head.  By 
successively  extracting  control  points  from  the  data,  we  were  able  to 
reproduce  and  display  the  head  quite  accurately. 


33 


-  Data  Reduction 

In  many  instances  it  is  possible  to  reduce  the  amount  of  data  needed  to 
properly  reconstruct  a  surface.  For  example,  consider  geographical  terrain. 
Terrain  that  is  relatively  flat  can  be  reconstructed  with  fewer  surface  patches 
that  can  mountainous  terrain.  The  problem  is  that  most  terrain  is  sampled 
at  discrete  intervals,  such  as  every  100  meters,  whether  it  is  flat  or  not.  By 
applying  some  form  of  an  Adaptive  Subdivision  Algorithm  [Ref.  5]  one  can 
reduce  the  amount  of  primary  and  secondary  storage  while  at  the  same  time 
provide  increased  performance  for  display. 

-  Lighting  Models 

Because  the  user  can  intercept  individual  polygons  comprising  the  surface,  it 
is  possible  to  subject  each  polygon  to  a  lighting  model.  While  we  have 
provided  a  simple  example  of  this,  more  sophisticated  lighting  models  could 
be  easily  integrated  through  these  parallel  functions. 

-  Realistic  3-D  Objects 

The  surfaces  of  many  vehicles  such  as  automobiles,  aircraft,  and  ships  can  be 
constructed  with  bicubic  surface  patches.  For  example,  constructing  an 
object  with  surface  patches  and  applying  a  lighting  and  shading  model,  one 
could  develop  an  ship  identification  training  system.  Such  a  training  system 
would  be  a  valuable  asset  in  military  training  environments,  allowing  the 
trainee  to  view  a  particular  class  of  ship  from  any  viewing  angle. 

2.  Improvement  of  Performance 


Real-time  computer  graphics  requires  efficient  algorithms  and  data 
structures.  While  these  functions  were  coded  to  be  as  efficient  as  possible,  while 
preserving  understandability,  there  is  always  room  for  improvement.  One 
suggestion  we  have  is  to  contact  the  developers  of  the  IRIS  graphics  package  for 
insights  into  improving  our  packages  performance.  Such  contact  may  provide 
access  to  low-level  graphic  system  routines  and  techniques  that  could  dramatically 
improve  performance. 
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B.  CONCLUSIONS 


This  study  introduced  the  reader  to  the  world  of  parametric  bicubic  surfaces. 
To  do  this,  we  provided  some  necessary  definitions,  terminology  and  mathematics. 
We  also  designed  and  implemented  a  set  of  software  functions  that  take 
advantage  of  the  information  and  given  them  to  the  reader  for  experimentation. 
The  benefit  that  can  be  derived  from  the  use  of  these  functions  can  only  be 
determined  by  the  passage  of  time. 
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APPENDIX  A 


FUNCTION  SPECIFICATIONS 


NAME 


ndefbasis  -  defines  a  basis  matrix 

SPECIFICATION 

ndefbasis(id,  mat) 
long  id: 

Matrix  mat; 

DESCRIPTION 

ndefbasis  allows  the  user  to  define  basis  matrices  for  use  in  the 
generation  of  patches,  matrix  is  saved  and  is  associated  with  id.  id 
may  then  be  used  in  subsequent  calls  to  npatchbasis. 


NAME 


npatchbasis  -  sets  current  basis  matrices 

SPECIFICATION 

npatchbasis(sid,  tid) 
long  sid,  tid; 

DESCRIPTION 

npatchbasis  sets  the  current  basis  matrices  (defined  by  ndefbasis  ) 
for  both  the  s  and  t  parametric  directions  of  a  surface  patch.  The 
current  s  and  t  bases  are  used  when  the  npatch  command  is  issued. 
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NAME 


npatchcurves  -  sets  number  of  curves  used  to  represent  a  patch 

SPECIFICATION 

npatchcurves  (scurves,  tcurves) 
long  scurves,  tcurves; 

DESCRIPTION 

npatchcurves  sets  the  current  number  of  s  and  t  curves  used  to 
represent  a  patch  as  a  wireframe. 


NAME 


npatchprecision  -  is  a  null  function. 

SPECIFICATION 

npatchprecision (ssegments,  tsegments) 
long  ssegments,  tsegments; 

DESCRIPTION 

npatchprecision  has  no  functionality  at  the  current  time.  It  is  used 
to  maintain  consistency  with  the  standard  IRIS  function 
patchprecision. 
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NAME 


npatch  -  draws  a  surface  patch 

SPECIFICATION 

npatch(geomx,  geomy,  geomz) 

Matrix  geomx,  geomy,  geomz; 

DESCRIPTION 

npatch  draws  a  surface  patch  using  the  current  npatchbasis  and 
npatchcurves.  The  shape  of  the  patch  is  determined  by  the  control 
points  specified  in  geomx,  geomy,  and  geomz. 


NAME 


Set  User  Routine  for  npatch  -  allows  the  user  to  specify  an 
intercept  function 

SPECIFICATION 

Set  User  Routine  for  npatch(fname) 
int  (*fname)(); 

DESCRIPTION 

Set  User  Routine  for  npatch  allows  the  user  to  set  up  a  function 
that  is  capable  of  intercepting  triangular  polygons  generated  during 
the  decomposition  of  a  surface  patch.  The  number  of  polygons 
generated  is  (  (seurves  -  1)  *  (tcurves  -  1)  *  2.  ). 
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NAME 


User  Routine  -  allows  the  user  to  turn  the  intercept  function  on  and 
off 

SPECIFICATION 

User  Routine(boolean) 
int  boolean; 

DESCRIPTION 

User  Routine  acts  like  a  switch  allowing  a  user-defined  intercept 
function  to  be  turned  on  and  off.  By  assigning  boolean  the  value  0  the 
intercept  function  is  turned  off.  Integer  values  other  than  0  cause  the 
intercept  function  to  be  turned  on. 


NAME 


Set  Default  Routine  for  npatch  -  resets  the  intercept  function 
to  a  system  defined  default 

SPECIFICATION 

Set  Default  Routine  for  npatchQ 

DESCRIPTION 

Set  Default  Routine  for  npatch  enables  the  user  to  choose  the 
system  defined  intercept  function  poly (3,  Triangle). 
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APPENDIX  B 


DEMONSTRATION  PROGRAMS 


/*  This  is  file  "basis. h"  */ 

^define  HERMITE  0 
#define  BEZIER  1 
#define  CARDINAL  2 
#define  BSPLINE  3 

/*  the  HERMITE  BASIS  MATRIX  */ 
Matrix  hermitematrix  =  { 

{  2.0,  -2.0,  1.0,  1.0  }, 

{  -3.0,  3.0,  -2.0.  -1.0  }, 

{  0.0.  0.0,  1.0,  0.0  }, 

{  1.0.  0.0.  0.0,  0.0  } 

}; 

/*  the  CARDINAL  BASIS  MATRIX  */ 
Matrix  cardinalmatrix  =  { 

{  -0.5,  1.5, -1.5,  0.5  }, 

{  1.0  ,  -2.5,  2.0,  -0.5  }, 

{ -0.5,  0.0,  0.5.  0.0  }, 

{  0.0.  1.0.  0.0,  0.0  } 

}; 

/*  the  BEZIER  BASIS  MATRIX  */ 
Matrix  beziermatrix  =  { 

{ -1.0,  3.0, -3.0,  1.0  }, 

{  3.0, -6.0,  3.0,  0.0  }, 

{ -3.0,  3.0,  0.0.  0.0  }, 

{  1.0.  0.0.  0.0,  0.0  } 

}; 

/*  the  B-SPLINE  BASIS  MATRIX  */ 
Matrix  bsplinematrix  =  { 

{  -1 .0/6.0,  3. 0/6.0,  -3. 0/6.0,  1.0/6.0  }, 
{  3. 0/6.0,  -6. 0/6.0,  3. 0/6.0,  0.0  }, 

{ -3.0/6.0,  0.0,  3.0/6.0,  0.0  }, 

{  1.0/6.0,  4.0/6.0,  1.0/6.0,  0.0  } 

}; 
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/*  This  is  file  "geom.h"  */ 

/*  set  up  the  geometry  matrix  of  x  coordinates  */ 

Coord  geomx[4][4]  =  { 

{  0.0.  100.0,  200.0.  300.0  }, 

{  0.0.  100.0.  200.0.  300.0  }. 

{  1000.0.  900.0,  800.0,  700.0  }, 

{  1000.0.  900.0,  800.0,  700.0  } 

}: 

/*  set  up  the  geometry  matrix  of  y  coordinates  */ 

Coord  geomy[4][4]  =  { 

{  400.0.  500.0,  600.0,  700.0  }, 

{  0.0.  200.0.  400.0,  600.0  }, 

{  0.0.  200.0.400.0,600.0  ), 

{  400.0.  500.0,  600.0,  700.0  } 

h 

/*  set  up  the  geometry  matrix  of  z  coordinates  */ 


Coord  geomz[4][ 

4]  =  { 

{ 

0.0. 

200.0. 

400.0, 

800.0 

{ 

0.0. 

200.0, 

400.0, 

800.0 

}, 

{ 

0.0. 

200.0, 

400.0, 

800.0 

}, 

{ 

0.0. 

200.0, 

400.0, 

800.0 

} 

}; 
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y*****************+********+****+***************************** 

Program  #1 

This  program  displays  two  wireframe  images  of  the  same 
surface  patch.  The  patch  drawn  in  the  color  YELLOW  is 
produced  by  the  standard  IRIS  patch  functions  and  the 
patch  drawn  in  the  color  RED  is  produced  by  the  parallel 
functions  we  have  developed. 

One  notices  that  the  patch  drawn  via  the  parallel  functions 
has  a  triangular  mesh  appearance  and  that  the  call  to  the 
npatchprecision  does  not  affect  the  displayed  image  as 
does  the  standard  IRIS  patchprecision  function.. 

+  *  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +++  +  +  +  +  +  +  +4 1  j 

^include  "gl.h"  /*  IRIS  graphics  library  */ 

^include  !,basis.hn 
^include  "geom-h" 

#define  S  CURVES  10 
#define  T  CURVES  10 

main() 

{ 

/*  Loop  variables  */ 
int  pi,  p2; 

/*  initialize  the  graphics  system  */ 
ginit() ; 

doublebuffer  () ; 
gconfig() ; 
cursoffQ; 

/*  set  up  the  viewing  parameters  */ 
ortho(0.0,  1023.0,  0.0,  1023.0,  -1023.0,  1023.0); 
viewport(0,  1023,  0,  767); 

/*  clear  the  graphics  screen  to  BLACK  */ 

color(BLACK); 

clear() ; 

/*  Associate  an  id  number  with  a  basis  matrix  */ 
defbasis(BEZIER,  beziermatrix) ; 
defbasis (CARDINAL,  cardinalmatrix) ; 
defbasis(BSPLINE,  bsplinematrix) ; 

ndefbasis(BEZIER.  beziermatrix); 
ndefbasis(C ARDINAL,  cardinalmatrix) ; 
ndefbasis(BSPLINE,  bsplinematrix) ; 
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/**  Specify  how  many  curves  in  each 
parametric  direction.  **/ 
patchcurves(S  CURVES  .  T  CURVES): 
npatchcurves(S  CURVES  .  T  CURVES); 

/**  Make  the  basis  matrices  different 
for  each  parametric  direction.  **/ 
patchbasis(BEZIER,  CARDINAL); 
npatchbasis(BEZIER,  CARDINAL); 

/**  Cycle  through  the  patch  changing  the 
precision  that  the  individual  curves 
comprising  the  patch  are  drawn.  **/ 
for(pl  =  10,  p2  =  100;  pi  <  100;  pi  +=  5,  p2  -=  5)  { 

/*  Draw  the  image  via  the  IRIS  functions  */ 

viewport(0.  511,  0,  767); 

color(BLACK); 

clear(); 

color(YELLOW); 
patchprecision(pl,  p2): 
patch(geomx,  geomy,  geomz); 

/*  Draw  the  image  via  the  parallel  functions  */ 

viewport (512,  1023,  0.  767); 

color(BLACK); 

clearQ; 

color(RED); 

npatchprecision(pl.  p2); 
npatch(geomx,  geomy,  geomz); 

/*  display  the  wireframe  images  */ 

swapbuffersQ: 

sleep(l): 

/*  clear  the  screen  */ 
color  (BLACK); 
clear(); 

} 

/*  clear  the  graphics  screen  and  exit  */ 

color  (BLACK); 

clearQ; 

gexitQ; 
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/*****************************  ******************************* 
Program  #2 

This  program  illustrates  how  the  triangles  formed  during 
a  surface  patch  decomposition  can  be  put  into  an  IRIS 
graphical  object  and  subsequently  displayed. 

# include  "gl.h"  /*  IRIS  graphics  library  */ 

^include  "basis. h" 

# include  "geom.h" 

#define  ON  1 

^define  S  CURVES  10 
#define  T  CURVES  10 

/*  Where  we  put  our  intercepted  triangles  */ 

Object  intercepted  object; 

main() 

{ 

/*  declare  the  intercept  function  */ 
int  intercept  function (); 

/*  initialize  the  graphics  system  */ 
ginit(); 

doublebuffer  (); 

gconfig(); 

cursoff(); 

/*  Make  the  intital  object.  */ 
makeobj(intercepted_object  =  genobj()); 
closeobj(); 

/*  set  up  the  viewing  parameters  */ 
ortho(0.0,  1023.0,  0.0,  1023.0,  -1023.0,  1023.0); 
viewport(0,  1023,  0,  767); 

/*  clear  the  graphics  screen  to  BLACK  */ 
color  (BLACK); 
clear  (); 

/*  Associate  an  id  number  with  a  basis  matrix  */ 
ndefbasis(BEZIER,  beziermatrix); 
ndefbasis (CARD INAL,  cardinalmatrix); 
ndefbasis(BSPLINE,  bsplinematrix); 

/**  Specify  how  many  curves  in  each 
parametric  direction.  **/ 
npatchcurves(S  CURVES  ,  T  CURVES); 
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/**  Make  the  basis  matrices  different 
for  each  parametric  direction.**/ 
npatchbasis(BEZIER,  BSPLINE); 

/**  Set  up  an  intercept  function  to  grab 
the  triangles  from  the  surface  patch 
and  turn  it  on.  **/ 

Set  User  Routine  for  npatch (intercept  function); 
User  Routine  (ON): 

/*  Call  the  npatch  function.  */ 
npatch  (geomx.  geomy,  geomz); 

/*  Display  the  surface  patch.  */ 
color(  YELLOW); 
callobj(intercepted  object); 
swapbuffersj); 
sleep(lG); 

/*  clear  the  graphics  screen  and  exit  */ 

color  (BLACK) : 

clear(); 

gexit(); 


} 


^tt***************************************************** 

This  is  the  user-defined  intercept  function  that 
handles  the  individual  triangles  generated 
during  the  decomposition  of  a  surface  patch. 

*******************************************************/ 

int  intercept  function(triangle) 
float  triangle[3]  [3]; 

{ 

/*  Open  up  the  object  and  put  in  the  triangles  */ 
editobj(interceptedobject); 

move(triangle[0][0],  triangle[0]  [l],  triangle[0]  [2]) ; 
draw(triangle[l]  [0],  triangle[l][l],  triangle[l]  [2]) ; 
draw(triangle[2]  [0],  triangle [2] [l] ,  triangle[2]  [2] ) ; 
draw(triangle[0]  [0] 9  trianglefo]  [l],  triangle[0]  [2] ) ; 
closeobjQ ; 

} 
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************************************************************* 
Program  #3 

This  programs  illustrates  how  the  intercept  function 
can  be  changed  as  the  program  runs. 


************************* 


************************************ 


/ 


^include  "gl.h"  /*  IRIS  graphics  library  */ 

^include  "basis.h" 

#include  "geom.h" 

/**  Define  a  new  type  that  is  a  pointer  to  a  function 
that  returns  an  integer.  **/ 
typedef  int  (*Function_pointer)  ()  ; 

^define  ON  1 

#define  S  CURVES  10 
^define  T  CURVES  10 

main() 

{ 

/*  Loop  variable  */ 
int  count; 

/*  Declare  an  array  of  pointers  to  functions  */ 
Function  pointer  intercept  functions[4] ; 

/*  declare  the  intercept  functions  */ 
int  intercept _functionl(); 
int  intercept _function2(); 
int  intercept_function3(); 
int  intercept  function4(); 

/**  initialize  the  array  of  pointers  to  the 
intercept  functions.  **/ 
intercept  functions[0]  =  intercept  functionl; 
intercept_functions[l]  =  intercept  function2; 
intercept_functions[2]  =  intercept  function3; 
intercept  functions [3]  =  intercept  function 4; 

/*  initialize  the  graphics  system  */ 
ginit(); 

doublebuffer  ( ) ; 
gconfig() ; 
cursoff(); 

/*  set  up  the  viewing  parameters  */ 
ortho(0.0,  1023.0,  0.0,  1023.0,  -1023.0,  1023.0); 
viewport(0,  1023,  0,  767); 
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/*  clear  the  graphics  screen  to  BLACK  */ 

color  (BLACK); 

clear(); 

/*  Associate  an  id  number  with  a  basis  matrix  */ 
ndefbasis(BEZIER,  beziermatrix); 
ndefbasis(  CARD  INAL,  cardinalmatrix); 
ndefbasis(BSPLINE.  bsplinematrix) ; 

/**  Specify  how  many  curves  in  each 
parametric  direction.  **/ 
npatchcurves(S  CURVES  ,  T  CURVES); 

/**  Make  the  basis  matrices  different 
for  each  parametric  direction.  **/ 
npatchbasis(BSPLINE,  CARDINAL); 

/*  Initially  use  the  default  intercept  function  */ 

Set  Default  Routine  for  npatch(); 

User  Routine  (ON); 

/*  Step  through  each  intercept  function  */ 
for(count  =  0;  count  <  4  ;  countH — h)  { 

/**  Set  up  an  intercept  function  to  grab 

the  triangles  from  the  surface  patch.  **/ 

Set  User  Routine  for  npatch (intercept  functions[count]); 

/**  Call  the  npatch  function  using  the  current 
intercept  function.  **/ 
npatch(geomx,  geomy,  geomz); 

/*  Display  what  the  intercept  function  did.  */ 

swapbuffers(); 

sleep  (2); 

/*  Clear  the  screen  and  do  another  one.  */ 

color  (BLACK) ; 

clear(); 

} 

/*  clear  the  graphics  screen  and  exit  */ 

color  (BLACK); 

clear(); 

gexitQ; 
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/********** **********************************  ************ 

This  intercept  function  draws  each  triangle  RED. 

************ **********************************♦****♦****/ 

int  intercept  functionl (triangle) 
float  triangle[3]  [3]; 

{ 

color(RED) ; 
poly (3,  triangle); 

} 


^* ******************************************************* 

This  intercept  function  draws  each  triangle  YELLOW. 

****** . *** . ***** . * . */ 

int  intercept _function2 (triangle) 
float  triangle[3] [3]; 

{ 

color  (YELLOW): 
poly(3,  triangle); 

} 


j  + ******************************************************* 

This  intercept  function  draws  each  triangle  GREEN. 

********************************************************  j 

int  intercept  function3 (triangle) 
float  triangle[3]  [3]; 

{ 

color(GREEN): 
poly(3,  triangle); 

} 


/+ ******************************************************* 

This  intercept  function  draws  each  triangle  BLUE. 

*************************************************+*****+/ 

int  intercept _function4(triangle) 
float  triangle[3] [3]; 

{ 

color  (BLUE); 
poly(3,  triangle); 

} 


48 


Program  #4 

This  program  illustrates  how  the  intercept  function 
can  be  combined  with  an  illumination  model  to 
provide  a  realistic,  illuminated,  three-dimensional 
solid-filled  curved  surface. 

+  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  f 


#include  "gl.h"  /*  IRIS  graphics  library  */ 

^include  ndevice.hfT 
#include  "basis. h,f 
#include  "geom.h" 

^define  ON  1 
#define  S  CURVES  25 
#define  T  CURVES  25 

main() 

{ 

/*  Declare  the  intercept  function.  */ 
int  light_poly(); 

/*  loop  variables  */ 
int  ij: 

/*  Initialize  the  graphics  system.  */ 
ginit(): 

singlebuffer  (); 
gconfig() ; 
cursoff(); 

/*  Clear  the  display.  * / 
color  (BLACK); 
clear(); 

/*  Set  up  new  viewing  parameters.  */ 
ortho(0.0,  1023.0.  0.0,  1023.0,  -1023.0,  1023.0); 
viewport(0,  1023,  0,  767); 

/*  Clear  drawing  area.  * / 

color(CYAN); 

clear(); 

/*  Set  for  hidden  surface  elimination.  */ 

setdepth(0x3FFF,  OxCOOO); 

zclear(); 

zbuffer(TRUE); 
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/*  Load  the  color  map  ramp  with  a  grey  scale.  */ 
for  (i  =  0;  i  <  256;  i-f — (-) 

{ 

mapcolor(8+i.  i,  i,  i); 

} 

/*  Associate  an  id  with  a  basis  matrix.  */ 
ndefbasis(BEZIER,  beziermatrix) : 
ndefbasis (CARDINAL,  cardinalmatrix); 
ndefbasis(BSPLINE,  bsplinematrix); 

/**  Provide  a  different  basis  for  each 
parametric  direction.  **/ 
npatchbasis(CARDINAL,  BEZIER); 

/**  Provide  the  number  of  curves  in 
each  parametric  direction.  **  f 
npatchcurves(S  CURVES  .  T  CURVES): 

/*  Set  up  the  intercept  function.  * / 

Set  User  Routine  for  npatch (light  poly); 

User  Routine(ON); 

while(TRUE) 

{ 

/*  Clear  the  z  buffer.  */ 
zclearQ; 

/*  Hold  display  if  MOUSE2  is  down.  */ 
if(getbutton(MOUSE2)) 

{ 

/*  Resume  when  MOUSE1  is  pressed  *•*/ 
while(!getbutton(MOUSEl))  ; 

} 

/*  Exit  when  MOUSEl  and  MOUSE2  and  MOUSE3  are  down.  */ 
if(getbutton(MOUSEl)  getbutton(MOUSE2)  &&  getbutton(MOUSE3)) 
break; 

/*  Clear  the  drawing  area.  */ 

color(CYAN); 

clear  (); 

/*  Set  the  current  color.  */ 
color  (BLACK); 

/*  Draw  initial  surface  patch.  */ 
npatch (geomx,  geomy,  geomz); 
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/**  Change  the  y  coordinates  to  get  a 
different  surface  patch  the  next 
time  we  draw  display.  **/ 
for(i=l;  i<3;  i++) 

{ 

for(j=0;  j<4;  j++) 

{ 

geomy[i][j]  =  (float)  (lrand48()  %  900); 

} 

} 

/*  Draw  the  surface  patch.  */ 
npatch(geomx,  geomy,  geomz): 

/**  Change  the  y  coordinate  values  to 
make  another  surface  patch.  **/ 
for(i=l;  i<3;  i++) 

{ 

for(j=0:  j<4;  j++) 

{ 

geomy[i][j]  =  (lrand48()%2)*geomy[i]  [j] : 

} 

} 

} 

/*  Clean  up  and  exit  the  program.  */ 
color  (BLACK) ; 
clearQ; 
gexit  () : 


/************************************************************ 

The  user-defined  intercept  function  used  to  grab  the 
triangles  generated  in  the  surface  patch  decomposition. 

int  light  poly(Triangle) 

Coord  Triangle[3]  [3]; 

{ 

/*  Put  each  triangle  through  an  illumination  model.  */ 

lightthepoly (Triangle,  3,  350.0,-1750.0,  350.0,  350.0,  1750.0,  350.0,  9,  264) 

} 
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************************************************************* 


/ 


lightpoly.c 

It  is  a  routine  that  computes  lighting  for  a  polygon  based 
upon  the  angle  between  the  Normal  vector  of  the  polygon 
and  the  direction  to  the  light  source. 

light  thepoly(xyz, ncoords, ax,ay,az, lx, ly,lz,colormin, colormax) 
xyz[][3]  =  floating  coords  of  the  polygon, 
ncoords  =  number  of  coordinates. 

ax,ay.az  =  interior  point  of  the  whole  object.  Used  to  determine 
outward  facing  normal  of  the  polygon.  This  is  the  same 
point  of  reference  that  would  be  used  for  backface 
polygon  removal. 

lx,ly.lz  =  vector  pointing  in  direction  of  the  light  source. 

colormin,  colormax  =  indices  used  for  the  colors  assigned  to  this 
polygon.  The  user  is  responsible  for  setting 
up  the  color  ramp. 

Note:  the  routine  also  puts  the  polygons  out  ordered  counterclockwise 
with  respect  to  the  interior  point  for  ease  of  backface  polygon 
removal. 

*************************************************************^ 
^include  <gl.h> 

#include  <math.h> 

^define  PIDIY2  1.570796327 
#define  CLOCKWISE  1 
#define  ROW  3 

light  thepoly  (xyz,  ncoords,  ax,  ay,  az,  lx,  ly,lz,colormin.  colormax) 

Coord  xyz[][3]; 

unsigned  int  ncoords; 

Coord  ax,ay,az:  /*  interior  point  of  the  whole  object.  */ 

Coord  lx,ly,lz;  /*  direction  to  the  light  source  */ 

int  col  or  min, col  or  max;  /*  color  min  /max  indices  */ 

{ 

/*  temp  coord  hold  */ 

Coord  *txyz; 

/*  loop  temps  */ 

register  unsigned  short  int  i  j; 
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/*  direction  test  function  */ 
int  npoly  orient (); 

/*  vectors  used  to  compute  the  polygon’s  normal  */ 

Coord  vl[3],v2[3]: 

/*  the  polygon’s  normal  */ 

Coord  normal[3]; 

/*  normal's  magnitude  */ 

Coord  normalmag; 

/*  light’s  magnitude  */ 

Coord  lightmag: 

/*  dot  product  of  N  and  L  */ 
double  dotprod; 

/*  angle  between  N  and  L  */ 
float  radians; 

/*  color  to  use  in  drawing  the  polygon  */ 
unsigned  short  int  colortouse; 

/*  allocate  memory  for  a  temporary  array  */ 

txyz  =  (Coord  *)  calloc  ((ncoords  *  3),  sizeof(Coord) ); 

/**  orient  the  polygon  so  that  its  counterclockwise  with  respect 
to  the  interior  point  **/ 

if(npoly  orient  (ncoords^xyz, ax, ay ,az)  ==  CLOCKWISE) 

{ 

/*  the  polygon  is  clockwise,  reverse  it.  */ 
for(i=0:  i  <  ncoords;  i=i  +  l) 

{ 

for(j=0;  j  <  ROW;  j=j+l) 

{ 

*(txyz  +  (i  *  ROW)  +  j)  =  xyz[ncoords-i-l] [j]; 

} 

} 

} 

else 

{ 

/*  no  need  to  reverse  */ 
for(i=0;  i  <  ncoords;  i=i+l) 

{ 

for(j=0;  j  <  ROW;j=j+l) 

{ 

*(txyz  +  (ROW  *  i)  +  j)  =  xyz[i][j]; 

} 

} 

} 
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/*  the  coordinates  are  ordered  counterclockwise  in  array  txyz  */ 

/**  compute  the  normal  vector  for  the  polygon  using  the  first 
three  vertices...  **/ 

/*  compute  the  first  vector  to  use  in  the  computation  */ 
vl[0]  =  *  (txyz  +  6)  -  *  (txyz  +  3);  /*  txyz[2][0]  -  txyz[l][0]  */ 

vl[l]  =  *  (txyz  +  7)  -  *  (txyz  +  4);  /*  txyz[2][l]  -  txyz[l][l]  */ 

vl[2]  =  *  (txyz  +  8)  -  *  (txyz  +  5);  /*  txyz[2][2]  -  txyz[l][2]  */ 

/*  compute  the  second  vector  to  use  in  computing  the  normal  */ 

v2[0]  =  *  (txyz  )  -  *  (txyz  +  3);  /*  txyz[0][0]  -  txyz[l][0]  */ 
v2[l]  =  *  (txyz  +  1)  -  *  (txyz  +  4);  /*  txyz[0][l]  -  txyz[l][l]  */ 
v2[2]  =  *  (txyz  +  2)  -  *(txyz  +  5);  /*  txyz[0][2]  -  txyz[l][2]  */ 

/*  the  normal  is  vl  x  v2  */ 
normaljO]  =  vl[l]*v2[2]  -  vl[2]*v2[l]; 
normal[l]  =  vl[2]*v2[0]  -  vl[0]*v2[2]: 
normal[2)  =  vl[0]*v2[l]  -  vl[l]*v2[0]: 

/*  compute  the  magnitude  of  the  normal  */ 
normalmag  =  sqrt(  (normal  [0]  *  normal  [0])  + 

(normal[l]*normal[l])  + 

(normal  [2]  *normal  [2] ) ) ; 

/*  compute  the  magnitude  of  the  light  */ 
lightmag  =  sqrt((lx  *  lx)  +  (ly  *  ly)  +  (lz  *  lz)); 

/*  compute  N  .  L  (normal  dot  product  with  the  light  source  direction)  */ 
dotprod  =  (double)  ( (normal [0]  *  lx) 

+  (normal[l]  *  ly) 

+  (normal[2]  *  lz)); 

/*  compute  the  unit  normal  */ 

dotprod  =  (double)  ((dotprod/ (normalmag  *  lightmag))); 

/*  dotprod  =  cos(theta)  of  the  angle  between  N  and  L. 

Convert  to  angle  in  radians  */ 
radians  =  acos(dotprod) ; 
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/*  compute  the  color  we  should  use  */ 
if(-PIDIV2  <=  radians  &&  radians  <=  PIDIV2) 

{ 

/*  if  the  angle  is  negative,  set  to  positive  */ 
if(radians  <  0.0) 

{ 

radians  =  -radians; 

} 

colortouse  =  ( (colormax-colormin)  /PIDIV2)  *  (PIDIV  2-radians)  +colormin; 

} 

else 

{ 

colortouse  =  colormin; 

} 

/*  set  the  color  */ 
color  (colortouse); 

/*  draw  the  poly  */ 
polf(ncoords,txyz): 

/*  free  up  memory  allocated  for  the  temporary  array  */ 
cfree(txyz,  (ncoords  *  3),  sizeof(Coord)); 


} 
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************************************************************ 


/ 


npoly  orient. c 

This  routine  determines  a  polygon’s  orientation 
with  respect  to  its  normal  and  a  reference  point. 
Orientation  is  either  clockwise  or  counter-clockwise. 
The  point  of  reference  must  not  lie  in  the  polygon’s 
plane. 


************************************************************ 


/ 


^include  <gl.h> 

^include  <math.h> 

int  npoly  orient(ncoords,xyz.xinside,yinside,zinside) 
unsigned  int  ncoords; 

Coord  xyz[][3]; 

Coord  xinside.  yinside,  zinside: 

{ 

/*  loop  temps  */ 

register  unsigned  short  int  i  j; 

/*  center  coordinate  of  the  polygon  */ 

Coord  center [3]; 

/**  vector  hold  locations  for  the  vectors  that  run 
from  the  center  coordinate  to  the  points  of  the 
polygon  **/ 

Coord  a[3],  b[3]; 

/**  points  on  line  containing  normal  that  are 
on  opposite  sides  of  the  plane  containing 
the  polygon.  **/ 

Coord  xn[3],  xmn[3j; 

/*  distance  to  point  n  from  pt  inside.  */ 
float  distton; 

/*  distance  to  point  -n  from  pt  inside.  */ 
float  disttomn; 

/*  the  normal  vector  computed  from  a  x  b  ’*'/ 

Coord  normal [3]; 

/*  compute  the  center  coordinate  of  the  polygon  */ 
center[0]  =  0.0; 
center[l]  =  0.0; 
center[2]  =  0.0; 
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for(i=0;  i  <  ncoords;  i-\ — b) 

{ 

for(j=0;  j  <  3;  j++) 

{ 

center[j]  -f=  xyz[i][j]; 

} 

} 

/*  divide  out  by  the  number  of  coordinates  */ 
for(j=0;  j  <  3;  j++) 

{ 

center[j]  =  center[j]/(float)ncoords; 

} 

/**  check  the  first  2  coordinates  of  the 
polygon  for  their  direction  **/ 

/**  compute  vector  a.  It  runs  from  the 
center  coordinate  to  coordinate  0  **/ 
for(j=0;  j  <  3:  j++) 

{ 

a  [j]  =  xyz[0][j]  -  center[j]; 

} 

/**  compute  vector  b.  It  runs  from  the 
center  coordinate  to  coordinate  1  **/ 
for(j=0;  j  <3;  j-\ — b) 

{ 

b[j]  =  xyz[l][j]  -  center  [j]; 

} 

/*  compute  a  x  b  to  get  the  normal  vector  */ 
normal[0]  =  a[l]’*‘b[2]  -  a[2]*b[l]; 
normalfl]  =  a[2]*b[0]  -  a[0]*b[2]; 
normal[2]  =  afoj^bfl]  -  a[l]*b[0]; 

/**  compute  point  n,  offset  pt  from  center  in 
direction  of  normal  **/ 
for(j=0;  j  <  3;  j++) 

{ 

xn[j]  =  centerfj]  +  normal[j]; 

} 

/**  compute  point  -n,  offset  pt  from  center 
in  opposite  direction  from  normal.  **/ 
for(j=0;  j  <  3;  j++) 

{ 

xmn[jj  =  center  [j]  -  normal[j]; 

} 
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/*  compute  the  distance  the  inside  pt  is  from  point  n  */ 
distton  =  sqrt((xn[0]  -  xinside)  *  (xn[0]  -  xinside)  + 

(xn[l]  -  yinside)  *  (xn[l]  -  yinside)  + 

(xn[2]  -  zinside)  *  (xn[2]  -  zinside)); 

/*  compute  the  distance  the  inside  pt  is  from  point  -n  */ 
disttomn  =  sqrt((xmn[Oj  -  xinside)  *  (xmn[Oj  -  xinside)  + 
(xmn[l]  -  yinside)  *  (xmn[l]  -  yinside)  + 

(xmn[2]  -  zinside)  *  (xmn[2]  -  zinside)); 

/**  if  the  dist(n)  <  dist(-n),  then  n  points  back  towards  the 
inside  point  and  is  on  the  same  side  of  the  plane  as  inside, 
a  x  b  is  then  clockwise.  **/ 
if(distton  <  disttomn) 

{ 

return(l);  /*  clockwise  */ 

} 

else 

{ 

return(O);  /*  counterclockwise  */ 

} 


} 
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APPENDIX  C  -  BENCHMARK  PROGRAM 


/************************************************************* 

This  is  the  BENCHMARK  PROGRAM  used  to  test  the 
relative  performance  of  the  standard  IRIS  functions 
and  the  parallel  routines  to  those  standard 
functions. 

********************************************************* ***4y 

^include  "gl.h"  /*  IRIS  graphics  library  */ 

^include  "basis. h,r 
^include  "geom.h" 

#define  IRIS  /*  Which  set  to  test  switch.  */ 

#define  MAX  TIMES  THRU  100 
#define  S  CURVES  25 

#define  T  CURVES  25 

fdefine  ONE  1 

main() 

{ 

int  times  thru; 

/*  initialize  the  graphics  system  */ 
ginit(): 

doublebuffer  () : 

gconfigQ; 

cursofif(); 

/*  clear  the  graphics  screen  */ 

color  (BLACK); 

clearQ; 

/*  set  up  the  viewing  parameters  */ 
ortho(0.0,  1023.0,  0.0,  1023.0,  -1023.0,  1023.0); 
viewport(0,  1023,  0,  767); 

/*  clear  the  graphics  screen  to  CYAN  */ 

color(CYAN); 

clear(); 
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#ifdef  IRIS 

/*  Use  the  standard  IRIS  functions  */ 
defbasis(BEZIER,  beziermatrix); 
defbasis(CARDINAL,  cardinalmatrix); 
defbasis(BSPLINE,  bsplinematrix); 

patchbasis(BEZIER,  BEZIER); 
patchcurves(S  CURVES  ,  T  CURVES); 
patchprecision(ONE,  ONE); 

#else 

/*  Use  the  parallel  functions  */ 
ndefbasis(BEZIER,  beziermatrix); 
ndefbasis(CARDINAL,  cardinalmatrix); 
ndefbasis(BSPLINE.  bsplinematrix); 

npatchbasis(BEZIER,  BEZIER); 
npatchcurves(S  CURVES  ,  T  CURVES); 
npatchprecision(ONE,  ONE); 
jjfendif 
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for(times_thru  =  0:  times  thru  <  MAX  TIMES  THRU:  times  thru-) — h)  { 

/*  clear  the  graphics  screen  to  CYAN  each  time  thru  */ 

color(CYAN); 

clear(); 

/*  draw  the  wireframe  surface  patch  in  BLACK  */ 
color(BLACK): 

#ifdef  IRIS 

/*  Using  the  IRIS  function  */ 
patch(geomx,  geomy,  geomz); 

#else 

/*  Using  the  parallel  function  */ 
npatch(geomx,  geomy,  geomz); 

fendif 

/*  display  the  wireframe  image  */ 
swapbuffers(); 

} 

/*  clear  the  graphics  screen  and  exit  */ 

color  (BLACK); 

clearQ; 

gexitQ; 


} 
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APPENDIX  D  -  PARALLEL  FUNCTIONS  SOURCE  CODE 


y*****+******+**+*++***+******************+*+**************** 

FILE  .  npatch.c 

Author .  Gary  W.  TAYLOR  (Captain  USMC) 

Date  .  1  December  1986 

Place  .  Naval  Postgraduate  School,  Monterey  CA 


Environment  . 

Silicon  Graphics.  Inc.,  IRIS  2400 
graphics  workstation,  UNIX  operating 
system  (GL2-W3.4). 


Purpose . 

The  following  C  source  code  provides  a  set 
of  "shadow"  routines  to  parallel  the  standard 
IRIS  2400  graphics  workstation  surface  patch 
routines.  These  parallel  routines  provide 
their  user  the  capability  to  generate 
solid-filled  parametric  bicubic  surface  patches. 


Notes  . 

As  of  the  current  date,  there  are  no  known 
side-effects  or  bugs  associated  with  using 
these  functions. 


Limitations  . 

These  functions  can  be  used  in  immediate  mode  only. 

+******++*++*++++*+**+++++++++++++++++****+*** 
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^include  "gl.h"  /*  IRIS  graphics  library  */ 

^include  "stdio.h" 

fdefine  OFF  0 

/*  Here  is  how  we  define  a  surface  point.  */ 
typedef  struct  { 

Coord  x; 

Coord  y; 

Coord  z; 

}  Point; 

/*  Structure  used  to  track  user  supplied  basis  matrices.  */ 
static  struct  list  elem  { 

/*  Matrix  id  number.  */ 

long  nid; 

/*  Pointer  to  the  basis  matrix.  */ 
float  *nmatrix; 

/*  Pointer  to  the  next  basis  matrix.  */ 
struct  list  elem  *next_elem  ptr; 


}; 
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/*  Used  in  forward  difference  computation  in  the  U  direction.  */ 
static  Matrix  Precision  Matrix  U; 

/*  Used  in  forward  difference  computation  in  the  V  direction.  */ 
static  Matrix  Precision  Matrix  V; 

/*  The  default  intercept  function.  */ 
static  int  Default  system  routine  (); 

/*  Pointer  to  the  currently  defined  U  basis  matrix.  */ 
static  float  *current  Ubasis; 

/*  Pointer  to  the  currently  defined  V  basis  matrix.  */ 
static  float  ^current  V  basis; 

/*  How  many  curves  in  the  U  direction.  */ 
static  long  UCURVES; 

/*  How  many  curves  in  the  V  direction.  */ 
static  long  VCURVES; 

/*  How  may  curve  segments  in  the  U  direction.  */ 
static  long  USEGMENTS; 

/*  How  may  curve  segments  in  the  V  direction.  */ 
static  long  VSEGMENTS; 

/*  Pointer  to  linked  list  of  user  supplied  basis  matrices.*/ 
static  struct  list  elem  *head_of_list  =  NULL; 

/*  Set  initial  user  routine  to  the  default.  */ 

static  int  (*  User  routine)  ()  =  Default  system  routine; 

/*  Initially  do  not  call  user's  function.  */ 
static  int  User  routine  is  on  —  OFF; 
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**************************************************************** 


Default  system  routine 

This  routine  is  what  the  system  automatically  does  if  the 
user  does  not  supply  a  particular  intercept  function. 


**************************************************************** 

static  int  Default  system  routine  (Triangle) 

Coord  Triangle[3] [3]; 

{ 

/*  Draw  the  polygon  with  a  standard  IRIS  function.  */ 
poly  (3,  Triangle); 

} 


/**************************************************************** 

Set  UserRoutinefornpatch 

This  function  allows  the  user  to  supply  an  intercept  function 
to  be  used  to  manipulate  the  triangles  that  are  generated  in 
the  decomposition  of  a  surface  patch. 

****************************************************************  j 

void  Set  User  Routine  for  npatch  (routine) 
int  (*routine)  (); 

{ 

/*  Save  the  pointer  to  the  user’s  function.  */ 

User  routine  =  routine; 

} 
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**************************************************************** 


/ 


User  Routine 

This  function  allows  the  user  to  turn  the  intercept  function 
on  or  off  at  will. 

**************************************************************** 

void  User  Routine  (boolean) 
int  boolean; 


/* 

If  boolean  =  0  then  the  intercept  function  is  turned  OFF. 
If  boolean  !=  0  then  the  intercept  function  is  turned  ON. 

v 


User  routine  is  on  =  boolean: 


} 


^/***************************************+************************ 


Set  Default  Routine  for  npatch 

This  function  allows  the  user  to  reset  the  intercept  routine 
to  the  same  routine  used  by  the  system. 


**************************************4 1************************* 

void  Set  Default  Routine  for  npatch  () 

{ 

/*  Point  to  the  default  intercept  function.  */ 

User  routine  =  Default  system  routine; 

} 
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**************************************************************** 


ndefbasis 

This  function  is  equivalent  to  the  IRIS  defbasis  function 
in  that  it  allows  the  user  to  define  a  basis  matrix  for 
use  in  the  generation  of  patches,  matrix  is  saved  and 
is  associated  with  id.  id  may  then  be  used  in  subsequent 
calls  to  npatchbasis. 


****************************************************************y 

void  ndefbasis  (id,  matrix) 
long  id: 

Matrix  matrix; 

{ 

/*  Special  processing  first  time  this  function  is  called.  */ 
static  int  has  been  called  =  FALSE; 

/*  Data  structure  pointer  for  new  entry.  */ 
struct  listelem  *new_elem_to_add; 

/*  Pointer  to  search  linked  list.  */ 
struct  listelem  *  walking  _ptr; 

/*  Pointer  to  a  copy  of  the  user  supplied  basis  matrix.  */ 
float  *pmatrix; 

/*  Loop  variables  */ 
int  row: 
int  column; 

/*  Get  memory  for  the  new  data  elements.  */ 

new  elem  to  add  =  (struct  list  elem  *)  malloc  (sizeof  (struct  list  elem)); 
pmatrix  =  (float  *)  calloc  (sizeof  (Matrix),  sizeof  (float)); 

/*  Make  a  copy  of  the  basis  matrix  passed  in  by  the  user.  */ 
for  (row  —  0;  row  <  4;  row+-f ) 
for  (column  =  0:  column  <  4;  column-f +)  { 

*  (pmatrix  +  (4  *  row)  +  column)  =  matrix[row]  [column]; 


} 


/*  Associate  the  user  supplied  id  to  this  basis  matrix.  */ 
new  elem  to  add  ->  nid  =  id; 

/*  Point  to  the  copied  basis  matrix.  */ 
new  elem  to  add  ->  nmatrix  —  pmatrix; 
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/**  Determine  how  to  add  this  information  into  the  linked  list 
of  basis  matrices.  **/ 

/*  Does  a  list  already  exist?  */ 
switch  (has  been  called)  { 

case  TRUE: 

/*  Point  to  the  beginning  of  the  list.  */ 
walking  ptr  =  head  of  list; 

/**  Traverse  the  list  looking  to  see  if  the  id  number 
already  exists.  **/ 

while  ((walkingptr  ->  nid  !=  id)  && 

(walking  ptr  ->  next  elem  ptr  !=  head  of  list))  { 

/*  Walk  through  the  linked  list.  */ 

walking  ptr  =  walking  ptr  ->  next  elem  ptr; 

} 


/*  id  already  exists  so  we  can  reuse  its  allocated  memory.  */ 
if  (walking  ptr  ->  nid  ==  id)  { 

/*  Get  rid  of  the  old  basis  matrix.  */ 

cfree  (walking  ptr  ->  nmatrix,  sizeof  (Matrix),  sizeof  (float)); 

/*  Point  to  the  replacement  matrix.  */ 
walking  ptr  ->  nmatrix  =  pmatrix; 

/*  Get  rid  of  the  un-needed  data  structure.  */ 
cfree  (new  elem  to  add,  1,  sizeof  (struct  list  elem)); 

} 

else  {  /*  The  id  does  not  exist  */ 

/**  Manipulate  the  pointers  to  add  the  new 
data  element  to  the  linked  list.  **/ 

new  elem  to  add  ->  next  elem  ptr  =  head  of  list  ->  next  elem  ptr; 
head  of  list  ->  next  elem  ptr  =  new  elem  to  add; 

} 

break; 


68 


case  FALSE: 


/*  No  linked  list  of  basis  matrices  exists  so  we  start  up  one.  */ 

/*  Create  the  pointer  to  the  front  of  the  list.  */ 
head  of  list  =  newelem  to  add; 
head  of  list  ->  next  elem  ptr  =  head  of  list: 

/*  Make  sure  we  do  not  do  this  again.  */ 
hasbeencalled  =  TRUE; 

break; 


default: 

break; 


} 


} 
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/+++++++*+++++++++++++++++*++++*++++*++++*+++++++++++++++++++++++ 

npatchbasis 

This  function  is  equivalent  to  the  IRIS  patchbasis  function 
in  that  it  sets  the  current  basis  matrices  for  both  the 
U  and  V  parametric  directions  of  a  surface  patch.  The  current 
U  and  V  bases  are  used  when  the  npatch  command  is  issued. 

++*+++++++*+*+*+*+++++**++*+*++++++++++++++++++*++++*+++++*+++++ I 

int  npatchbasis  (uid,  vid) 
long  uid,  vid; 

{ 

struct  list  elem  *walking_ptr; 

/*  ERROR:  no  linked  list  of  basis  matrices.  */ 
if  (head  of  list  ==  NULL)  { 

fprintf  (stderr,  "Opatchbasis:  no  basis  matrices  definedO); 
exit  (-1); 


} 

/*  Traverse  the  list  looking  for  the  desired  U  basis  matrix.  */ 
walking  ptr  =  head  of  list; 
while  ((walking  ptr  ->  nid  !=  uid)  && 

(walking  ptr  ->  next  elem  ptr  !=  head  of  list))  { 

walkingptr  =  walking  ptr  ->  next  elem  ptr; 


} 

if  ((walking  ptr  ->  nid  !=  uid)  && 

(walking  ptr  ->  next  elem  ptr  ==  head  of  list))  { 

/*  ERROR:  U  basis  matrix  does  not  exist  in  the  linked  list.  * / 

fprintf  (stderr,  ,fOpatchbasis:  undefined  U  basis  matrix  %d0?  uid); 
exit  (-1); 


} 

current  U  basis  =  walking  ptr  ->  nmatrix; 
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/*  Traverse  the  list  looking  for  the  desired  V  basis  matrix.  */ 
walking  ptr  =  head  of  list; 
while  ((walking  ptr  ->  nid  !=  vid) 

(walking  ptr  ->  next  elem  ptr  !=  head  of  list))  { 

walking  ptr  =  walking  ptr  ->  nextelemptr; 


} 

if  ((walking  ptr  ->  nid  !=  vid)  &&: 

(walking  ptr  ->  next  elem  ptr  ==  head  of  list))  { 

/*  ERROR:  V  basis  matrix  does  not  exist  in  the  linked  list.  */ 

fprintf  (stderr,  "Opatchbasis:  undefined  V  basis  matrix  %d0.  vid); 
exit  (-1)  : 


} 

current  V  basis  =  walking  ptr  ->  nmatrix; 


} 
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t*************************************************************** 


/ 


npatchcurves 

This  function  is  similar  to  the  IRIS  patchcurves  command, 
ucurves  and  vcurves  set  the  subdivision  parameters  used  in 
decomposing  the  surface  patch.  An  individual  patch  will  be 
decomposed  into  a  (ucurves  -1)  *  (vcurves  -1)  grid  with 
each  grid  generating  two  triangular  polygons. 

void  npatchcurves  (ucurves,  vcurves) 
long  ucurves,  vcurves; 

{ 

/*  Prevent  an  inappropriate  number  of  curves.  */ 

UCURVES  =  (ucurves  <  2  ?  2  :  ucurves  -  1); 

VCURVES  =  (vcurves  <  2  ?  2  :  vcurves  -  1); 

/**  Set  up  the  Precision_Matrix  U  used  in  the  forward  difference 
along  the  U  direction.  **/ 

Precision  Matrix  U[0][0]  =  6.0  /  (float)  (UCURVES  *  UCURVES  *  UCURVES); 

Precision  Matrix  U[1][0]  =  6.0  /  (float)  (UCURVES  *  UCURVES  *  UCURVES); 

Precision  Matrix  U[l][l]  =  2.0  /  (float)  (UCURVES  *  UCURVES); 

Precision  Matrix JJ [2] [0]  =  1.0  /  (float)  (UCURVES  *  UCURVES  *  UCURVES); 

Precision  Matrix  U[2][l]  =  1.0  /  (float)  (UCURVES  *  UCURVES); 

Precision  Matrix_U[2]  [2]  =  1.0  /  (float)  (UCURVES): 

Precision  Matrix  U[3][3]  =  1.0; 

/**  Set  up  the  Precision  Matrix  V  used  in  the  forward  difference 
along  the  V  direction.  **/ 

Precision  Matrix  V[0][0]  =  6.0  /  (float)  (VCURVES  *  VCURVES  *  VCURVES); 

Precision  Matrix  V[1][0]  =  6.0  /  (float)  (VCURVES  *  VCURVES  *  VCURVES); 

Precision  Matrix  V(l][l]  =  2.0  /  (float)  (VCURVES  *  VCURVES); 

Precision  Matrix  V [2] [0]  =  1.0  /  (float)  (VCURVES  *  VCURVES  *  VCURVES); 

Precision  Matrix  V[2][l]  =  1.0  /  (float)  (VCURVES  *  VCURVES); 

Precision_Matrix_V[2] [2]  =  1.0  /  (float)  (VCURVES); 

Precision_Matrix_V[3][3]  =  1.0; 


} 
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**************************************************************** 


npatchprecision 

This  function  is  similar  to  the  IRIS  patchprecision  routine 
but  is  only  used  to  maintain  consistency  with  the  IRIS 
routines.  Its  results  are  not  used  by  any  other  function 
in  the  suite. 


**************************************************************** 

void  npatchprecision  (usegments.  vsegments) 
long  usegments.  vsegments; 

{ 

/*  Prevent  an  inappropriate  number  of  segments.  */ 

USEGMENTS  =  (usegments  <  2  ?  2  :  usegments  -  1); 

VSEGMENTS  =  (vsegments  <  2  ?  2  :  vsegments  -  1); 

} 
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/++++*+++++++++++++++++++++++*+++++++++++++++++++++++++++++++++++ 

nspeckle 

This  function  finishes  computation  and  hands  off  each  triangle 
to  an  intercept  function. 

**********************************************+*****************^ 

static  void  nspeckle  (Coord  array) 

Point  *  Coord  array [4]; 

{ 

register  Point  *  Patch  array; 

/*  Get  enough  memory  to  hold  all  the  points  that  will  be  generated.  */ 
Patch  array  =  (Point  *)  calloc  ((UCURVES  +  1)  *  (YCURVES  +  1), 
sizeof  (Point)); 

{ 

/*  Pointer  to  a  particular  point.  */ 
register  Point  *  Where; 

register  unsigned  int  totalpoints; 
register  unsigned  int  point  count; 
register  unsigned  int  t  count; 
register  unsigned  int  count; 
register  unsigned  int  Row; 
register  unsigned  int  Column; 

/*  Used  in  generating  points  on  the  surface.  */ 

Matrix  control  matrix; 

/*  Intermediate  matrix  to  hold  mathematical  results.  */ 

Matrix  interl; 

/*  For  every  point  in  the  the  U  parametric  direction.  */ 
for  (pointcount  =  0,  total  points  =  0; 

pointcount  <=  UCURVES;  point  countH — |-)  { 

/*  Build  a  control  matrix  for  the  current  curve.  */ 
for  (count  =  0;  count  <  4;  count++)  { 

Where  =  (Point  *)  (Coord  array [count]  +  point_count); 

c on trol_matrix[ count]  [0]  =  Where  ->  x; 
control_matrix[count]  [l]  =  Where  ->  y; 
control_matrix[count][2]  =  Where  ->  z; 

} 
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/**  Generate  the  matrix  to  compute  the  forward  difference  on. 

The  forward  difference  matrix  is  equal  to 

Precision  Matrix  V  x  current  V  basis  x  control  matrix.  **/ 

pushmatrix  (); 
loadmatrix  (controlmatrix) ; 
multmatrix  (current  V  basis): 
multmatrix  (Precision  Matrix  Y); 
getmatrix  (interl); 
popmatrix  (); 

/*  Generate  the  points  on  the  curve  in  the  V  direction  */ 
for  (t  count  =  0;  t  count  <=  VCURVES:  t  count-| — K  total  points++)  { 

(Patch  array  +  total  points)  ->  x  =  interl[3] [0] ; 

(Patch  array  +  total  points)  ->  y  =  interl [3]  [l] ; 

(Patch  array  +  total  points)  ->  z  =  interl[3] [2]; 

/*  Do  the  forward  difference.  */ 
for  (Row  =  3;  Row  >  0;  Row—) 
for  (Column  =  0:  Column  <  4;  Column++) 
interl  [Row]  [Column]  =  interl  [Row]  [Column]  +  interl[Row  -  l][Column]; 


} 


} 


} 
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{ 


/*  Place  to  put  a  triangle  to  send  out  to  user  */ 

Coord  Triangle_l[4]  [3]; 

register  Point  *  Where: 
register  unsigned  int  Row; 
register  unsigned  int  Column; 

/*  Decompose  the  patch  into  its  individual  triangles.  */ 
for  (Row  =  0;  Row  <  UCURVES;  Row++) 
for  (Column  =  0;  Column  <  VCURVES;  Column++)  { 

Where  =  (Patcharray  +  ((VCURVES  +  1)  *  Row)  +  Column); 

Triangle  l[0][0]  =  Where ->  x; 

Triangle  l[0][l]  =  Where  ->  y; 

Triangle  1[0]  [2]  =  Where  ->  z; 

Triangle  l[l][0]  =  (Where  +  1)  ->  x; 

Triangle_l[l]  [l]  =  (Where  +  1)  ->  y; 

Triangle  l[l][2]  =  (Where  +  1)  ->  z; 

Triangle_l[2][0]  =  (Where  +  VCURVES  +  1)  ->  x; 

Triangle  1  [2] [1]  =  (Where  +  VCURVES  +  1)  ->  y; 

Triangle  1[2]  [2]  =  (Where  +  VCURVES  +  1)  ->  z; 

Triangle_l[3] [0]  =  (Where  +  VCURVES  +  2)  ->  x; 
Triangle“l[3][l]  =  (Where  +  VCURVES  +  2)  ->  y; 

Triangle  1  [3]  [2]  =  (Where  +  VCURVES  +  2)  ->  z; 

/*  Does  the  user  have  an  intercept  routine?  */ 
if  (  User  routine  is  on)  { 

/*  Yes  */ 

(*_User_routine)  (&Triangle_l[0] [0]); 

(*  User  routine)  (&Triangle_l[l] [0]); 

} 

else  { 

/*  No  user  routine,  so  we  use  the  default.  */ 

Default  system  routine  (&Triangle_l  [0]  [0] ) ; 

Default  system  routine  ( &Triangle_l [lj  [0] ) ; 

} 

} 


/*  Return  the  memory  used  to  the  system.  */ 

cfree  (Patch  array,  (UCURVES  +  1)  *  (VCURVES  +  1),  sizeof  (Point)); 

} 
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^4444444444444444444444444444444444444444444444444444444444444444 

npatch 

This  function  rearranges  the  input  matrices  into  a  form  readily 
used  in  computing  points  along  the  four  curves  defined  by  those 
matrices.  It  then  computes  points  for  each  curve  in  the  U 
direction  using  the  technique  of  forwards  differencing  of  a 
matrix.  Using  these  points  we  can  then  generate  points  along 
the  a  curve  in  the  V  direction. 

4444444444444444444444444444444444444444444444444444444444444444^ 

void  npatch  (geomx,  geomy,  geomz) 

Coord  geomx[4][4],  geomy [4] [4],  geomz[4][4]; 

{ 

register  Point  *  Coord_array[4]; 

/*  One  control  matrix  for  each  curve.  */ 

Matrix  ctrl  ptsl; 

Matrix  ctrl_pts2; 

Matrix  ctrl_pts3; 

Matrix  ctrl  pts4; 
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/*  Load  the  appropriate  control  matrix  for  each  curve.  */ 


/*  Curve  1  */ 

ctrl_ptsl[0]  [0]  =  geomx[0][0]; 
ctrl_ptsl[0]jl]  =  geomyjojjoj; 
ctrl_ptsl[0]  [2]  =  geomz[0][0]; 

Ctrl  ptsl[l] [0]  =  geomx[l][0]; 
ctrl_ptsl[l][l]  =  geomyjl]  jo]; 
ctrl_ptsl[l]  [2]  =  geomz[l][0]; 

ctrl_ptsl[2]  [0]  =  geomx[2][0]; 
ctrl_ptsl[2]  [l]  =  geomy[2][0]; 
ctrl_ptsl[2]  [2]  =  geomz[2][0]; 

ctrl_ptsl[3]  [0]  =  geomx[3][0]; 
Ctrl  ptsl[3]  [l]  =  geomy[3]  [0]; 
ctrl_ptsl[3]  [2]  =  geomz[3][0]; 

/*  Curve  2  */ 

ctrl_pts2[0]  [0]  =  geomx[0][l]; 
ctrl_pts2[0]  [1]  =  geomyjojjlj; 
ctrl_pts2[0]  [2]  =  geomz[0][l]; 

ctrl_pts2[l][0]  =  geomx[l][l]; 
ctrl_pts2[l][l]  =  geomyjl]  jlj; 
ctrl_pts2[l][2]  =  geomz[l]  [1]; 

Ctrl  pts2[2][0]  =  geomx[2][l]; 
ctrl_pts2[2]  [1]  =  geomy|2][lj; 
ctrl_pts2[2]  [2]  =  geomz[2]  [l]; 

Ctrl  pts2[3]  [0]  =  geomx[3][l]; 
ctrl_pts2[3]  jl]  =  geomy[3][lj; 
ctrl_pts2[3]  [2]  =  geomz[3][l]; 
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/*  Curve  3  */ 

Ctrl  pts3[0][0]  =  geomx[0]  [2] ; 
Ctrl  pts3[0][l]  =  geomy[0] [2] ; 
ctrl_pts3[0]  [2]  =  geomz[0]  [2] : 

Ctrl  pts3[l][0]  =  geomx[l][2]; 
ctrl_pts3[l][l]  =  geomy[l][2]; 
ctrl_pts3[l]  [2]  =  geomz[l][2]; 

Ctrl  pts3[2]  [0]  =  geomx[2][2]; 
Ctrl  pts3[2][l]  =  geomy[2][2]; 
Ctrl  pts3[2][2]  =  geomz[2][2]; 

ctrl_pts3[3]  [0]  =  geomx[3][2]; 
ctrl_pts3[3][l]  =  geomy[3][2]; 
ctrl_pts3[3][2]  =  geomz[3]  [2]: 

/*  Curve  4  */ 

Ctrl  pts4[0][0]  =  geomx[0] [3] ; 
Ctrl  pts4[0][l]  =  geomy[0]  [3] ; 
ctrl_pts4[0][2]  =  geomz[0]  [3] : 

ctrl_pts4[l][0]  =  geomx[l]  [3] ; 
Ctrl  j)ts4[l]  [l]  =  geomyjl]  [3]; 
ctrl_pts4[l][2]  =  geomz[l][3]; 

Ctrl  pts4[2]  [0]  =  geomx[2]  [3]; 
ctrl_pts4[2][l]  =  geomy[2][3]; 
ctrl_pts4[2]  [2]  =  geomz[2] [3]; 

Ctrl  pts4[3][0]  =  geomx[3][3]; 
ctrl_pts4[3][l]  =  geomy[3][3]; 
ctrl_pts4[3][2]  =  geomz[3][3]; 
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{ 


register  Point  *Where: 

register  unsigned  int  Row; 
register  unsigned  int  Column; 
register  unsigned  int  point  count; 
register  unsigned  int  count; 

/*  An  array  of  pointers  to  our  control  matrices.  */ 
float  *matrix  pointer[4]; 

/*  Matrix  used  to  hold  mathematical  results.  */ 

Matrix  interl; 

/*  Initialize  the  array  */ 
matrix  pointer[0]  =  (float  *)  Ctrl  ptsl; 
matrix  pointerfl]  =  (float  *)  ctrl  pts2; 
matrix  pointer[2]  =  (float  *)  ctrl_pts3; 
matrix  pointer[3]  =  (float  *)  ctrl  pts4; 

/*  Get  enough  memory  to  hold  the  points.  */ 

Coord  arrayfO]  =  (Point  *)  calloc  (UCURVES  +  1.  sizeof  (Point)); 

Coord  arrayfl]  =  (Point  *)  calloc  (UCURVES  +  1,  sizeof  (Point)); 

Coord  array[2]  =  (Point  *)  calloc  (UCURVES  +  1,  sizeof  (Point)): 

Coord  array[3]  =  (Point  *)  calloc  (UCURVES  +  1,  sizeof  (Point)); 

/*  For  each  curve.  * / 
for  (count  =  0;  count  <  4;  countH — f-)  { 

/*  Generate  the  matrix  used  in  the  forward  difference  for  this  curve.  */ 
pushmatrix  (); 

loadmatrix  (matrix  pointerfcount]); 
multmatrix  (current _U  basis); 
multmatrix  (Precision  Matrix  U); 
getmatrix  (interl); 
popmatrix  (); 
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/*  For  each  curve  generate  points  in  the  U  parametric  direction.  */ 
for  (point  count  =  0;  point  count  <=  UCURVES;  point  count++)  { 

Where  =  (Point  *)  (Coord  array[count]  +  point  count) ; 

Where  ->  x  =  interl[3] [0] ; 

Where  ->  y  =  interl [3] [lj; 

Where  ->  z  =  interl[3]  [2]; 

/*  Do  the  forward  difference.  */ 
for  (Row  =  3;  Row  >  0;  Row—) 
for  (Column  =  0;  Column  <  4;  Column  +  +  ) 
interl  [Row]  [Column]  =  interl  [Row]  [Column]  +  interl  [Row  -  l]  [Column]; 


} 


} 


} 

/*  Call  function  to  finish  computations  and  display.  */ 
nspeckle  (Coord  array); 

/*  Return  the  memory  used  back  to  the  system.  */ 
cfree  (Coord  array [0],  (UCURVES  +  1),  sizeof  (Point)); 
cfree  (Coord _array[l],  (UCURVES  +  1),  sizeof  (Point)); 
cfree  (Coord_array[2],  (UCURVES  +  1),  sizeof  (Point)); 
cfree  (Coord  array [3],  (UCURVES  +  l),  sizeof  (Point)); 


} 
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